Introductory Tutorials

This section provides an overview of the basic functionalities offered by CRISP. Each tutorial introduces a specific feature of the package, guiding you through typical use cases.

First verify the installation by importing CRISP and its dependencies:

import CRISP
print(f"CRISP version: {CRISP.__version__}")

CRISP is organized into two main subpackages:

• simulation_utility: Tools for preparing and processing simulation data

• data_analysis: Methods for analyzing simulation results

The following tutorials are organized according to these subpackages.

Simulation Utility

This section covers utilities for manipulating and preparing simulation data.

Atomic Indices

CRISP allows for the classification of atoms based on custom-defined indices, making it easier to analyze specific subsets within your simulation data.

1. Classify Atom Indices

   from CRISP.simulation_utility.atomic_indices import run_atom_indices
   
   # Path to your ASE trajectory file
   file_path = "./wrapped_traj.traj"
   
   # Output folder to save the indices
   output_folder = './indices_new/'
   
   # Define the custom cutoffs dictionary
   custom_cutoffs = {
       ("O", "H"): 1.2,
       ("Si", "O"): 1.8,
       ("Al", "Si"): 3.2,
       ("O", "O"): 3.0
   }
   
   # Run the atom_indices function and save the results with a specific frame index
   run_atom_indices(file_path, output_folder, frame_index=10, custom_cutoffs=custom_cutoffs)
   

   Output:

   Analyzing frame with index 10 (out of 21 frames)
   Length of O indices: 432
   Length of H indices: 120
   Length of Si indices: 168
   Length of Al indices: 24
   Outputs saved.
   Saved cutoff indices for O-H to ./indices_new/cutoff/O-H_cutoff.csv
   Saved cutoff indices for Si-O to ./indices_new/cutoff/Si-O_cutoff.csv
   Saved cutoff indices for Al-Si to ./indices_new/cutoff/Al-Si_cutoff.csv
   Saved cutoff indices for O-O to ./indices_new/cutoff/O-O_cutoff.csv
   

Atomic Trajectory Visualization

Visualize atomic trajectories in 3D to understand motion and structural changes.

from CRISP.simulation_utility.atomic_traj_linemap import plot_atomic_trajectory

# Path to trajectory file
traj_file = "./wrapped_traj.traj"

# Output directory for the visualization
output_dir = "./atomic_traj_linemap/o593_atom_trajectory.html"

# Select specific atom for trajectory analysis
selected_atoms = [593]

# Generate interactive 3D visualization of trajectories
fig = plot_atomic_trajectory(
    traj_path=traj_file,
    indices_path=selected_atoms,
    output_dir=output_dir,
    frame_skip=1,
    plot_title="Oxygen Atom Movement Analysis",
    atom_size_scale=1.2,
    show_plot=True
)

Output:

Loading trajectory from ./wrapped_traj.traj (using every 1th frame)...
Loaded 21 frames from trajectory
Using 1 directly provided atom indices for trajectory plotting
Simulation box dimensions: [24.34499931 24.34499931 24.34499931] Å
Analyzing atom types in first frame (total atoms: 744, max index: 743)...
Found 4 atom types: Si, Al, O, H
Plot has been saved to ./atomic_traj_linemap/o593_atom_trajectory.html/trajectory_plot.html

Visualisation Output:

This interactive 3D visualization allows you to:

• Rotate, zoom, and pan to explore atomic trajectories

• Track the movement of selected atoms over time

• Visualize the trajectory path with customizable atom sizes

• Understand atomic motion patterns within the simulation box

Subsampling

Extract representative structures from a trajectory using Farthest Point Sampling.

from CRISP.simulation_utility.subsampling import subsample

# Extract representative structures from trajectory
all_frames = subsample(
    traj_path="./wrapped_traj.traj",
    n_samples=10,
    index_type="all",
    file_format="traj",
    frame_skip=1,
    plot_subsample=True,
    output_dir="./Subsmapling"
)

print(f"Selected {len(all_frames)} representative structures")

Output:

Saved convergence plot to ./Subsmapling/subsampled_convergence.png
Saved 10 subsampled structures to ./Subsmapling/subsample_wrapped_traj.traj

Visualisation Output:

The convergence plot shows the distance between each sampled structure and its nearest neighbor, illustrating how the algorithm selects maximally diverse structures from the trajectory.

Error Analysis

Perform statistical error analysis on time-correlated simulation data using different methods.

Example 1: Position Data Analysis

from CRISP.simulation_utility.error_analysis import autocorrelation_analysis
import numpy as np

# Load position data
data_positions = np.load("./error/positions.npy")

# Analyze using autocorrelation method
res_positions = autocorrelation_analysis(
    data_positions,
    plot_acf=True,
    max_lag=500
)

print(res_positions)

Output:

{'mean': array([11.84336219,  6.56230374,  6.34512439]),
 'acf_err': array([0.11042688, 0.0483816 , 0.06882431]),
 'std': array([0.21002227, 0.09201757, 0.13089782]),
 'tau_int': 69.11286151958006,
 'optimal_lag': 109}

Visualization Output:

Example 2: Energy Data Analysis

from CRISP.simulation_utility.error_analysis import autocorrelation_analysis, block_analysis
import numpy as np

# Load energy data from log file
data_energy = np.loadtxt("./error/md_20k.log", skiprows=1, usecols=2)

# Analyze using autocorrelation method
acf_error = autocorrelation_analysis(data_energy, plot_acf=True)

# Analyze using blocking method
block_error = block_analysis(data_energy, convergence_tol=0.001, plot_blocks=False)

print(acf_error)
print(block_error)

Output:

{'mean': -3065.5796212000005,
 'acf_err': 0.0054549762052233325,
 'std': 0.6834300195116669,
 'tau_int': 0.318542675191111,
 'optimal_lag': 9}
{'mean': -3065.5796212000005,
 'block_err': 0.02208793396091249,
 'std': 0.6834300195116669,
 'converged_blocks': 32}

Visualization Output:

Interatomic Distance Calculation

Calculate and save distance matrices between atoms for further analysis.

from CRISP.simulation_utility.interatomic_distances import distance_calculation, save_distance_matrices

# Path to trajectory file
traj_path = "./wrapped_traj.traj"
frame_skip = 10
index_type = ["O"]  # Focus on oxygen atoms

# Calculate full and subset distance matrices
full_dms, sub_dms = distance_calculation(traj_path, frame_skip, index_type)

# Save the calculated distance matrices
save_distance_matrices(full_dms, sub_dms, index_type, output_dir="distance_calculations_zeo")

Output:

Distance matrices saved in 'distance_calculations_zeo/distance_matrices.pkl'

This utility calculates distance matrices between atoms, accounting for periodic boundary conditions, and saves the results for later use in clustering or other analyses.

Data Analysis

This section covers methods for analyzing simulation results and extracting physical insights.

Contact and Coordination Analysis

CRISP provides tools to analyze both coordination environments and dynamic contacts between atoms.

1. Coordination Analysis

from CRISP.data_analysis.contact_coordination import coordination

# Path to trajectory file
filename = "./wrapped_traj.traj"

# Define target and bonded atoms for coordination analysis
target_atoms = "O"
bonded_atoms = ["O"]

# Define custom cutoffs for coordination
custom_cutoffs = {('Al', 'O'): 3.5, ('O', 'O'): 2.5}

# Perform coordination analysis
cn = coordination(filename,
                 target_atoms, bonded_atoms, custom_cutoffs,
                 skip_frames=10, plot_cn=True, output_dir="./CN_data")

Output:

Interactive coordination distribution chart saved to ./CN_data/CN_distribution.html

Visualization Output:

The visualizations show both: - The distribution of coordination numbers across all analyzed frames, revealing dominant coordination environments - The time evolution of coordination numbers throughout the trajectory, showing structural changes over time

The visualization shows the distribution of coordination numbers for the selected atom types, allowing you to identify dominant coordination environments and their frequencies throughout the simulation trajectory.

2. Contact Analysis

from CRISP.data_analysis.contact_coordination import contacts

# Path to trajectory file
filename = "./wrapped_traj.traj"

# Define target and bonded atoms for contact analysis
target_atoms = "O"
bonded_atoms = ["O"]

# Define custom cutoffs for contacts
custom_cutoffs = {('Al', 'O'): 3.5, ('O', 'O'): 2.5}

# Perform contact analysis
sub_dm, cal_contacts = contacts(
    filename, target_atoms, bonded_atoms, custom_cutoffs,
    skip_frames=1,
    plot_distance_matrix=True,
    plot_contacts=True,
    time_step=50.0*1000,  # fs
    output_dir="./Contacts_data")

Output:

Interactive contact heatmap saved to ./Contacts_data/O_heatmap_contacts.html
Interactive contact analysis chart saved to ./Contacts_data/average_contact_analysis.html
Static contact analysis chart saved to ./Contacts_data/average_contact_analysis.png
Interactive distance heatmap saved to ./Contacts_data/O_heatmap_distance.html

Visualization Output:

The contact analysis provides multiple visualizations:

• A summary of average contact statistics and their distribution

• A distance heatmap showing the spatial relationships between atoms

• A time heatmap showing the persistence of contacts throughout the trajectory

These visualizations enable the identification of persistent contacts and transient interactions throughout the simulation trajectory.

Hydrogen-Bonding Analysis

Analyze hydrogen bond networks and dynamics in your simulation.

from CRISP.data_analysis.h_bond import hydrogen_bonds

# Perform hydrogen bond analysis
h_bonds_both_plots = hydrogen_bonds(
    traj_path="./wrapped_traj.traj",
    frame_skip=1,
    acceptor_atoms=["O"],
    angle_cutoff=120,
    mic=True,
    output_dir="./H_Bond_Data",
    time_step=50*1000,
    plot_count=True,
    plot_heatmap=True,
    plot_graph_frame=True,        # Generate frame-specific plot
    plot_graph_average=True,      # Generate average plot
    graph_frame_index=10          # Use frame 10 instead of default 0
)

Output:

Hydrogen bond count plot saved to './H_Bond_Data/h_bond_count.png'
H-bond structure 2D histogram saved to './H_Bond_Data/h_bond_structure.png'
Generated and saved 196 unique donor/acceptor atom indices to ./H_Bond_Data/donor_acceptor_indices.npy
Generating hydrogen bond network visualizations for frame 10...
Interactive correlation matrix saved as './H_Bond_Data/hbond_correlation_matrix_frame_10.html'
Figure saved as './H_Bond_Data/hbond_network_frame_10.html'
Generating average hydrogen bond network visualization...
Interactive correlation matrix saved as './H_Bond_Data/hbond_correlation_matrix_average.html'
Figure saved as './H_Bond_Data/hbond_network_average.html'

Visualization Output:

The hydrogen bond analysis provides multiple complementary visualizations:

• Time series of hydrogen bond counts throughout the trajectory

• Heatmap showing hydrogen bond structure and persistence over time

• Interactive network graph showing hydrogen bonding patterns for a specific frame (frame 10)

• Correlation matrix visualizing the strength and frequency of hydrogen bonds between specific atoms

These visualizations work together to provide a comprehensive understanding of hydrogen bonding networks, enabling both qualitative pattern recognition through the network graph and quantitative analysis through the correlation matrix. The analysis identified 196 unique donor/acceptor atom pairs and provides both frame-specific and time-averaged hydrogen bonding information.

Radial Distribution Function (RDF)

Perform Radial Distribution Function analysis to investigate spatial relationships between atoms.

from CRISP.data_analysis.prdf import analyze_rdf

# Path to trajectory file
traj_file = "./wrapped_traj.traj"

# RDF parameters
rmax = 10.0         # Maximum radius
nbins = 50          # Number of bins for histogram
frame_skip = 1      # Analyze every frame
output_dir = "custom_ase"
output_filename = None  # Auto-generate filename
use_prdf = False    # Calculate total RDF (not partial)
atomic_indices = None

# Perform RDF analysis
data_rdf = analyze_rdf(
    use_prdf=use_prdf,
    rmax=rmax,
    traj_path=traj_file,
    nbins=nbins,
    frame_skip=frame_skip,
    output_filename=output_filename,
    atomic_indices=atomic_indices,
    output_dir=output_dir,
    plot_prdf=True  # Generate plots
)

Output:

Data saved in 'custom_ase/rdf_total.pkl'
Static plot saved in 'custom_ase/rdf_total_plot.png'
Interactive animation saved to 'custom_ase/rdf_total_animation.html'
GIF animation saved to 'custom_ase/rdf_total_animation.gif'
Animation saved in 'custom_ase/rdf_total_animation.gif'

Visualization Output:

Once Loop Reflect

The RDF analysis provides both static and interactive visualizations:

• A static plot showing the time-averaged radial distribution function

• An interactive animation displaying the evolution of the RDF across multiple frames

These visualizations help identify characteristic distances between atoms and structural features such as coordination shells and their changes throughout the simulation.

Mean-Square Displacement (MSD)

Calculate and visualize Mean-Square Displacement to analyze diffusion in two steps.

Step 1: Calculate MSD values from trajectory

from CRISP.data_analysis.msd import calculate_save_msd

# Path to trajectory and indices files
traj_file = "./SiAl15/nvt.traj"
indices_file = "./SiAl15/indices_needed/ex_fram_ox.npy"
timestep = 50 * 10  # fs

# Calculate MSD values and save to CSV
msd_values, msd_times = calculate_save_msd(
    traj_path=traj_file,
    timestep_fs=timestep,
    indices_path=indices_file,
    output_file="msd_results.csv",
    frame_skip=10
)

Output:

Loaded 2200 frames after applying frame_skip=10
Loaded 72 atom indices
Using timestep: 500 fs
Calculating MSD using parallel processing...
MSD data has been saved to traj_csv_detailed/msd_results.csv

Step 2: Analyze MSD data and calculate diffusion coefficient

from CRISP.data_analysis.msd import analyze_from_csv
import pandas as pd

# Load the MSD data CSV file
df = pd.read_csv("./traj_csv_detailed/msd_results.csv")
print(f"Total data points in file: {len(df)}")

# Calculate diffusion coefficient from MSD data
D, error = analyze_from_csv(
    csv_file="./traj_csv_detailed/msd_results.csv",
    fit_start=0,
    fit_end=len(df),
    with_intercept=True,
    plot_msd=True,
    plot_diffusion=True,
    n_blocks=5
)

Output:

Total data points in file: 2200
Loaded MSD data from ./traj_csv_detailed/msd_results.csv
Diffusion Coefficient: 4.77e-07 cm²/s
Error in the Fit: 3.80e-09 cm²/s
Diffusion coefficient evolution plot saved to: ./traj_csv_detailed/msd_results_3D_diffusion_evolution.png
MSD Analysis Results (Standard Fit):
Diffusion Coefficient: D = 4.7665e-07 ± 3.8004e-09 cm²/s (0.8%)

Visualization Output:

The MSD analysis provides comprehensive quantitative information about particle diffusion:

• Diffusion coefficient (D): 4.77e-07 cm²/s with statistical error (3.80e-09 cm²/s)

• Relative error: 0.8%, indicating high precision in the measurement

• Visual representation: MSD vs. time with linear fit showing the diffusive behavior

• Evolution plot: 3D visualization showing how the diffusion coefficient converges over different time windows

This two-step approach allows for efficient analysis of large trajectories by first extracting the MSD data and then performing detailed statistical analysis without reprocessing the trajectory. The block analysis (n_blocks=5) provides additional confidence in the diffusion coefficient measurement by examining its convergence behavior.

Clustering Analysis

Identify atomic clusters using advanced clustering algorithms. CRISP offers both single-frame and trajectory-based clustering analysis.

1. Single-Frame Clustering

Analyze clusters in a specific frame of your trajectory:

from CRISP.data_analysis.clustering import analyze_frame
import numpy as np
import os

# Path to trajectory and indices files
traj_file = "./SiAl15/nvt.traj"
indices_file = "./SiAl15/indices_needed/ex_fram_ox.npy"

# Clustering parameters
threshold = 3.0
min_samples = 3
output_dir = "SiAl15_Frame_clustering"

os.makedirs(output_dir, exist_ok=True)

# Load atom indices
atom_indices = np.load(indices_file)

# Create analyzer instance and perform clustering
analyzer = analyze_frame(
    traj_path=traj_file,
    atom_indices=atom_indices,
    threshold=threshold,
    min_samples=min_samples,
    metric='precomputed',
    custom_frame_index=-1  # Analyze last frame
)

# Perform clustering analysis
results = analyzer.analyze_structure(output_dir=output_dir)

Output:

Saving results to directory: SiAl15_Frame_clustering
3D visualization saved to SiAl15_Frame_clustering/nvt_clusters.html

Number of Clusters: 10
Number of Outliers: 9
Silhouette Score: 0.2932
Average Cluster Size: 6.30
Cluster Information:
  Cluster 0: 14 points
  Cluster 1: 9 points
  Cluster 2: 4 points
  Cluster 3: 7 points
  Cluster 4: 10 points
  Cluster 5: 5 points
  Cluster 6: 3 points
  Cluster 7: 4 points
  Cluster 8: 4 points
  Cluster 9: 3 points
Detailed frame data saved to: SiAl15_Frame_clustering/frame_data.txt
Full analysis data saved to: SiAl15_Frame_clustering/single_frame_analysis.pkl

Visualization Output:

The 3D visualization shows the spatial distribution of clusters in the selected frame, with each cluster represented by a different color and unclustered atoms shown separately.

2. Trajectory-Based Clustering

Analyze clusters throughout a trajectory to observe their evolution:

from CRISP.data_analysis.clustering import analyze_trajectory, save_analysis_results, plot_analysis_results
import os
import numpy as np

# Path to trajectory and indices files
traj_file = "./SiAl15/nvt.traj"
indices_file = "./SiAl15/indices_needed/ex_fram_ox.npy"

# Clustering parameters
threshold = 3.0
min_samples = 3
skip_frames = 10  # Analyze every 10th frame
output_dir = "SiAl15_Traj_clustering"
output_prefix = "SiAl15_traj_clusters"

os.makedirs(output_dir, exist_ok=True)

# Analyze trajectory
analysis_results = analyze_trajectory(
    traj_path=traj_file,
    indices_path=indices_file,
    threshold=threshold,
    min_samples=min_samples,
    frame_skip=skip_frames,
    output_dir=output_dir,
    save_html_visualizations=True  # Save HTML visualizations of first and last frames
)

# Save and plot results
pickle_file = save_analysis_results(
    analysis_results=analysis_results,
    output_dir=output_dir,
    output_prefix=output_prefix
)

plot_analysis_results(pickle_file, output_dir=output_dir)

Output:

Loaded 72 atoms for clustering from ./SiAl15/indices_needed/ex_fram_ox.npy
Loading trajectory from ./SiAl15/nvt.traj (using every 100th frame)...
Loaded 220 frames from trajectory
Analyzing 220 frames...
3D visualization saved to SiAl15_Traj_clustering/first_frame_clusters.html
3D visualization saved to SiAl15_Traj_clustering/last_frame_clusters.html
Trajectory analysis complete: 220 frames processed
Analysis results saved to directory: SiAl15_Traj_clustering
Analysis plot saved to: SiAl15_Traj_clustering/SiAl15_traj_clusters_plot.png

Visualization Output:

The trajectory-based clustering analysis provides comprehensive insights into cluster evolution:

• Temporal Evolution: Analysis of 220 frames (every 100th frame) showing how clusters change over time

• 72 Atoms Analyzed: Focused analysis on oxygen atoms from the framework

• Comparative Visualizations: Interactive 3D plots of first and last frames showing structural changes

• Statistical Summary: Evolution plot tracking cluster statistics throughout the trajectory

This analysis is particularly valuable for understanding:

• Long-term structural stability and reorganization

• Phase transitions and nucleation events

• Self-assembly processes in molecular dynamics simulations

• Framework flexibility and dynamic behavior in zeolite-like materials

The combination of temporal analysis with 3D visualizations provides both quantitative metrics and qualitative insights into the clustering behavior of atoms throughout the simulation.

Volumetric Density

Create 3D volumetric density maps to visualize the spatial distribution of atoms throughout a trajectory.

from CRISP.data_analysis.volumetric_atomic_density import create_density_map
import os

# Path to trajectory and indices files
traj_file = "./SiAl15/nvt.traj"
indices_file = "./SiAl15/indices_needed/ex_fram_ox.npy"

projection_options = [True, False]

for show_projections in projection_options:
    output_dir = "./SiAl15_Density"
    os.makedirs(output_dir, exist_ok=True)

    proj_text = "with_projections" if show_projections else "no_projections"
    output_file = f"SiAl15_density_{proj_text}.html"

    create_density_map(
        traj_path=traj_file,
        indices_path=indices_file,
        frame_skip=100,
        threshold=0.0,
        opacity=0.8,
        absolute_threshold=False,
        save_density=True,
        omit_static_indices=indices_file,
        show_projections=show_projections,
        projection_opacity=0.8,
        save_projection_images=show_projections,
        output_dir=output_dir,
        output_file=output_file
    )

Output:

Loading trajectory from ./SiAl15/nvt.traj (using every 100th frame)...
Loaded 220 frames, 72 selected indices
Extracting selected atom positions from trajectory...
Creating density grid...
Density data saved to: ./SiAl15_Density/nvt_density_data.npz
Creating visualization with relative threshold=0.0, opacity=0.8
Density range: 0.0 to 29.0 counts
Adding 2D projections of density data...
XY projection saved to: ./SiAl15_Density/SiAl15_density_with_projections_xy_projection.png
YZ projection saved to: ./SiAl15_Density/SiAl15_density_with_projections_yz_projection.png
XZ projection saved to: ./SiAl15_Density/SiAl15_density_with_projections_xz_projection.png
Skipping atom visualization since projections are enabled
Visualization saved as HTML file: ./SiAl15_Density/SiAl15_density_with_projections.html
Loading trajectory from ./SiAl15/nvt.traj (using every 100th frame)...
Loaded 220 frames, 72 selected indices
Extracting selected atom positions from trajectory...
Creating density grid...
Density data saved to: ./SiAl15_Density/nvt_density_data.npz
Creating visualization with relative threshold=0.0, opacity=0.8
Density range: 0.0 to 29.0 counts
Omitting 72 custom indices from static structure visualization
Visualization saved as HTML file: ./SiAl15_Density/SiAl15_density_no_projections.html

Visualization Output:

The volumetric density analysis provides comprehensive visualization of atomic distribution:

• 3D Density Maps: Interactive volumetric visualizations showing where oxygen atoms are most frequently located

• 220 Frames Analyzed: Statistical sampling using every 100th frame for robust density calculations

• 72 Oxygen Atoms: Focused analysis on extra-framework oxygen atoms

• Dual Visualization Modes: - With projections: 2D cross-sectional views (XY, YZ, XZ) for detailed analysis - Without projections: Clean 3D view with static framework atoms visible

This analysis is particularly valuable for understanding:

• Identifying preferred locations of mobile species in zeolite frameworks

• Understanding diffusion pathways and bottlenecks

• Visualizing guest molecule distribution patterns

• Analyzing adsorption sites and their accessibility

Example Jupyter Notebooks

For detailed examples and interactive tutorials, refer to the Jupyter notebooks included in the package:

• Basic Usage Examples: https://github.com/Indranil17/CRISP_HOST/blob/main/example/CRISP_Introdcutory_Tutorial.ipynb

These notebook provide step-by-step examples with visualizations to help you understand how to use CRISP effectively.